Boundary Conditions 


To create a stellar model, one must solve four interrelated differ- 
ential equations, which describe how the radius, pressure, luminos- 
ity, and temperature change with mass. To do this, we need four 
boundary conditions (to start off each differential). The boundary 
conditions for mass and luminosity are obvious: 


rin =0)=0 (4.1.1) 


L(m = 0) =0 (4.1.2) 


Unfortunately, we do not have any constraints for the central tem- 
perature and pressure. For these variables, we must apply boundary 
conditions at the surface of the star. To first order, we will define 
the stellar “surface” as the place where the optical depth, 7 = 1. At 
this location, m = Mr, and we can apply the boundary conditions 





P(m = Mr) = 0 (4.1.3) 


T(m= Mr) =0 (4.1.4) 


This approximation isn’t too bad, considering the temperatures and 
pressures in the interior of the star. A slightly better approximation 
is to use the equation 


Lr = 40 R oT (4.1.5) 





(which actually defines effective temperature) and realize that the 
light we see must come from some finite optical depth 


ia kpdr=R | pdr 
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and note that 


P(R) =| apdr = go | pdr 
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Setting these two equations equal yields 
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where, again, T ~ 1. (Actually, it’s 2/3.) 





(4.1.6) 


An even better approximation can be obtained using relations be- 
tween surface temperature, pressure, luminosity, and radius that 
come from the calculations of stellar atmospheres. Grids of model 
stellar atmospheres exist that give values for P and T at their base, 
as a function for log g and £. By interpolating in these grids, bound- 





ary conditions for P and T can be provided. 


Expansions for the Center 


A quick inspection of the equations of stellar structure demonstrates 
that they cannot be evaluated at M = 0, due to powers of r and 
M in the denominator. For these locations, you need to use an 
expansion of each equation. For the continuity equation (2.1.3), 
this means 


1 1 
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If we assume that very close to the center, p ~ p®, and we apply the 
boundary condition r = 0 at M = 0, then 
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or 





r= (2 (4.3.1) 
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Similarly, if the gradient in the energy sources is small, (2.3.3) ex- 
pands to 


dL = (en —e&v +e) dM => L=(&— epte) M (4.3.2) 
with its boundary condition. 


For the momentum equation (2.2.9), we can assume that the center 
of the star is in mechanical equilibrium, so 
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where we have substituted for radius using (4.3.1). Integrating this 








yields 
Papos T [ M-3dM 
or 
Papias (A w (4.3.3) 


In the case of radiative energy transport, we can expand the tem- 
perature equation (3.1.5) to 
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where we have again substituted using (4.3.1) and (4.3.2). Keeping 
everything but the mass constant, we then get 
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yielding 
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Finally, for convective energy transport, we re-write (2.4.4) as 
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which implies 
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Solving the Equations 


The equations of stellar structure are not analytic; they must be 
solved numerically. Furthermore, because the boundary conditions 
for the equations aren’t all at the center, the differentials can’t just 
be integrated outward. 





The key to computing a stellar model is to divide a star up into K 
mass shells, with 10 < K < 1000. The distribution of these shells 
is somewhat of an art form: in energy generating zones, where the 
temperature and density must be known with high accuracy, and 
near the surface, where the gradients are steep, the mass fraction 
must be small, ~ 1074Mr. Zones as large as 0.1M- may be used 
elsewhere. In general, static stars need fewer shells than rapidly 
evolving stars. 


THE MATCHING POINT METHOD 


The original method of creating a numerical stellar model is the 
Matching Point, or Fitting method. In this procedure, one starts 
with the two center boundary conditions (r = £ = 0), but guesses 
the central temperature and pressure. ‘The star is then integrated 
outward. At the same time an inward integration is performed, 
starting with the known surface boundary conditions (P, and T;) 
and guesses about R and £r. One then attempts to match the 
integrations somewhere in the middle of the star. 





In general, of course, the two integrations won’t match up: the 
guesses for the surface radius and luminosity and the central tem- 
perature and pressure will be wrong. At the matching point, the 
equations for the inward and outward integrations will differ by an 
amount Y. However, suppose the initial model is not too far off. 


Let x; be the initial guesses for R, Lr, Te, and P.. (The index 
i, of course, goes from 1 to 4.) Perturb one of the initial guesses, 
such that x; = x; + dx;. Although the model still won’t agree at 
the matching point, you can compare the old values of Y; with the 
new values Y/. In other words, you obtain an estimate for 2% 
for each equation. Now compute similar partial derivatives for the 
other variables. You can now estimate the perturbations you need 
to match the inward and outward integrations, i.e., to make Y; = 0. 


This is done via 4 independent linear equations of the form 


OY, OY; OY OY; 
—— | 0 — ]ô —— ]óð — | dv, = -Y; 
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By solving the equations for the perturbations, ôx;, you can estimate 


the magnitude of the true perturbations needed to match things up. 
After a few iterations, the model should converge. 





NUMERICAL TRICKS 


Because of the large range of numbers involved, and the possibility of 
numerical roundoff in computers, there are a number of useful tricks 
that are sometimes used (and/or needed) in model calculations. 





e Instead of using mass as the independent variable, sometimes 
one defines a new variable €, such that M = M(¿) and0<€ <1. 
This changes the equations of stellar structure slightly, i.e., mass 
conservations becomes 
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but makes the program more flexible, in that by changing M(6), 
you can immediately adjust to different mass distributions within a 
star. (For example, the distribution of mass within a main sequence 
star is completely different from that of a giant star.) By changing 
to €, the calculations can be made more linear. 


e In computing derivatives by finite differences, some round off 
error may occur, especially if the variables are changing by many 
orders of magnitude. To minimize these errors, variables are gener- 
ally made as linear as possible. For pressure and density, this means 
changing variables again, t.e., substituting 





p= pi/4 q = p'i/? To= LJE 


Again, these substitutions change the character, but not the sub- 
stance of the equations. 


e In computing numerical derivatives, the most straightforward 
method is to just take the arithmetic average of quantities computed 
on each side of a shell. However, if the quantity being measured is 
increasing exponentially, an arithmetic average may not be appro- 
priate: a geometrical average may be better. Unfortunately, there 


is no one prescription for this problem, other than to use variables 
that are as linear as possible and to keep the shells thin. 


e When computing the perturbations ôr, dP, ôL, and ôT, it is 
often times useful to “sneak up” on a solution, and limit your max- 
imum correction to some small value. For example, if the equations 
predict that ôr should be increased by 5%, it may be better to cor- 
rect the value by only ~ 1%. This will ensure that the perturbations 
stay in the linear regime. Without such a limit, the non-linearities in 
the equations may cause you to overcorrect, and continually bounce 
back and forth over the true solution. 








The Henyey Method 


The trouble with using the Matching Point method is that each 
model must be computed from scratch. While this may work for 
calculating a single model (say, a single main sequence star), it is 
inefficient for computing a series of models, as is done when mod- 
eling stellar evolution. Thus, the method of choice for computing a 
series of stellar models is the Henyey method, after Henyey, Forbes, 
& Gould 1964, Ap. J., 139, 306. (Other names for it include the 
relaxation method, the difference method, the shell method, and 
the generalized Newton-Raphson method.) To use the method, one 
first must start with a guess as to the structure of the star. (Usu- 
ally, this is done by first computing a main-sequence model via the 
Matching Point Method.) Hopefully, this guess is not too far off — 
if it is, then the model will not coverge. 











To see how the Henyey method works, we can start by re-writing all 
the spatial differentials as finite differences. In other words, (2.1.3), 
(2.2.9), (2.3.3), and (2.4.4) are written in the form of 
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The equations of stellar structure then become 
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V= A;j4 (4.4.4) 


where the superscript j refers to the number of the shell, and A; j = 
0. Note that all the variables to the left of the equal signs are 
functions of M, t, r, P, £L, and T, and they are (usually) evaluated 
half-way between shell j and j + 1. Thus (4.4.1) is 
pith _ pl 1 1 1 
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Note also that the variables p, cp, 6, and Vaq are known through 
the equation of state, €n, €v, and « are functions of p, T, and atomic 
physics, and the time derivatives can be written in terms of known 
values from the previous model, t.e., 


dPi Pi — P’ 

a 4.4. 
dt At a 
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(If the second derivative term is needed, it is written via the velocity, 
which adds another variable and equation to the mix. The time 
steps for mechanical calculations, however, must be short, so that 
no shell has Ar/At exceeding the sound speed.) 


Next, we have the four special equations for the center of the star. 
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(4.4.11) 
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where again, the values A;,; = 0. Note that in deriving (4.4.8) and 
(4.4.9) we have applied two boundary conditions, r(0) = £(0) = 0. 


Finally, we have the two outer boundary conditions: 
P*® = Bp (4.4.13) 


cS Br (4.4.14) 





where Bp and Br are either zero, or obtained from L* and r* 
using a look-up table of stellar atmospheres. (Recall, however, that 


normally the interior integration isn’t carried out all the way to 
the surface; instead a match point is chosen interior to the region 
where partial ionization, super-adiabaticity, and changing opacities 
become important.) 


Using the above equations, and the initial guess of the results, we 
can proceed with calculating a stellar model. Suppose our initial 
guess of the stellar model is wrong. In that case, the values A; ; in 
equations (4.4.1) through (4.4.4) and (4.4.8) through (4.4.13) will 
not be zero. We will therefore need to change our values of r, P, L, 
and T by or, ôP, ôL, and ôT. By how much? If we are not too far 
off, then linear theory should hold, and, for example 
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There are 4 equations like this for each mass shell. Except for shell 
1, each equation contains 8 partial derivatives. (The innermost con- 
tains only six partials, since the values of r and £ at the center are 
known; see (4.4.8) through (4.4.12)). In addition, we have the two 
equations (with four partials each) that come from the outer bound- 
ary condition. Now note: each partial derivative can be computed 
directly using the values of the initial guess, and the values of A; j 
are already known. Thus, the only unknowns in (4.4.15) are the 
values of ôr, ôP, ôL, and ôT. 








We can write the series of (4.4.15)-type equations in matrix form. 
Suppose we have a (minimal) three shell model. The relations be- 
tween the partial derivatives and A; j, Bp, and Br is 
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where H is called the Henyey matrix. The unknown vector contain- 
ing the correction values is then 


U=H'A (4.4.17) 





In fact, since the Henyey matrix is rather sparse (with many zeros, 
except along the diagonal), it’s fairly easy to find efficient algorithms 
to compute its inverse. Once U is known, these correction values 
can be used to improve the original guess, and the procedure can 
be repeated. After several iterations, the values A; j, Bp, and Br 
may approach zero, giving you the structure of the star. 


Digression: Forward and Back Differencing Techniques 


There are two ways to numerically solve a system of differential 
equations. The first method is through “forward differencing,” or 
“explicit” numerical integration. With this technique, values for in- 
tegrating step n + 1 are estimated directly from the known values 
of the previous step. In the alternative “back differencing” or “im- 
plicit” integration method, the values at step n + 1 are estimated 


using the n + 1 values themselves, via numerical inversion. 





For example, consider the simple differential equation 
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which has the analytical solution 
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In the forward differencing method, this equation is 
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which translates to 
Xni1 = aAt+ X, (1 — bAt) (4.2.1) 


while in the back differencing method, the equation is 


Aagi E A 
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Now note that (4.2.1) and (4.2.2) behave very differently as the time 
step becomes large, i.e., At — oo. With forward differencing, the 
numerical solution diverges; with back differencing, X,141 — a/b, 
which is the analytical result. This is characteristic of the two ap- 
proaches. In general, forward differencing solutions are unstable, 
and can only be used with very small time steps (2.e., for the above 
example, bAt < 1). Back differencing methods are much more reli- 
able, and should be used whenever possible. 





Digression: Solving Initial Value Problems 


A “standard” way to solve a set of differential equations with known 
initial values is through fourth-order Runge-Kutta integration. To 
appreciate the method, let us solve the simple differential equation 
dy 
ni 
using a Taylor expansion. If we expand y(x) about x) = 0, the 
expression for the function becomes 





y(x ) = y(æo) +y (eo) As + E Aa ro) Ag? 2 Y (20) an (A1.2) 


where Ax = (x — xo). From the boundary condition, we know that 
y(xo) = 1, hence from (A1.1), 


gro) = f (0, Yo) = 20 tyo =O+1=1 (A1.3) 


Next, we take the derivative of both sides of (A1.1), to get 
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From (A1.3) we know that y’(xo) = 1, hence (A1.4) gives us 
iy (eo) =1+1=2 (A1.5) 


Similarly, if we take the next derivative of the function, 
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and evaluate it at x = xp using (A1.5), we get y” (xo) = 2. This 
process can be extended as far out as needed. In the present case, 
the evaluation of y(x) yields 


y(x) = 1+ Ax + (Ax)? + (Ar) PO- (A1.6) 


where O ~ Az is the error term. 


Note here, that had we stopped the expansion after the first term 
(i.e., keeping the relation linear), the expression for y(x) would have 
been 


y(x) = y(xo) + y' (xo) Az = y(z0) + f (20, yo)Ax+O*%  (A1.7) 


This is just simple forward-differencing solution (or Euler’s method). 
Obviously, from (A1.7), we can expect the errors associated with the 
technique to be ~ (Az)? per step. If the entire interval is divided 
into 1/Az steps, then we expect the total error in the numerical 
integration to be ~ Ag. 








It is easy to see why Euler’s method is not very accurate — the 
slope used to determine the new value of y always comes from the 
beginning of the interval, and thus is always wrong. A better esti- 
mate would come from using the slope of the function at the middle 
of the interval, i.e., 


y(2) = ylzo) + H 


Of course this function cannot be evaluated, since we cannot eval- 
uate y'(x) without knowing y(x). However, we can estimate y(x) 
from Euler’s method, and then use the predicted value of y’(x) to 
produce a “corrected” value of y(x). In other words 
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y(x) = (zo) + es 


(This procedure can be repeated again if necessary, although after a 
while, it’s better just to cut down the step size.) The accuracy of this 
modified Euler technique is ~ (Azx)?. To see this, simply compare 
(A1.8) with the Taylor series expansion (A1.2), while substituting 
a forward-differencing expression for y” (xo), i.e., 
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As above this error is per step; the total error of the integration is 
Ax?. 


The next improvement comes by combining the results of several 
estimates along the interval and weighting each one differently. For 
instance, suppose we generalize (A1.8) to 


y(x) = y(ao) + akı + bk2 where 
kı = y' (z0, y0) Az ko = y (£o + aAz, yo + Bky)Ax (A1.11) 


where a and b are the weights of the two estimates, œ describes 
where on the interval the estimates are made, and 8 is a correction 


factor for the estimate made using kı. When written out fully with 
the substitution f(x,y) = y’, (A1.11) becomes 
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If we expand the last term as a two-dimensional Taylor series about 
T0, Yo, then 
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This expression can then be compared to the Taylor expansion for 
y(x) given in (A1.2) 
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Thus, a+b = 1, ab = 1/2, and 6b = 1/2. These are the coefficients 
for the second-order Runge-Kutta method; note that since there are 
three equations and four unknowns, more than one set of coefficients 
will work. 





The fourth-order Runge-Kutta method works the same way, except 
that four intermediate estimates are made. The equations for this 
method are 
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In general, the method gives results with accuracies of the order 
of (Ax) per step, and (Az)* over the entire integration range. 
Because of its accuracy, and relatively few computational steps, it 
is perhaps the most widely used technique for solving initial-value 
differential equations. 





A final improvement, given by Press in Numerical Recipes, is to 
use fourth-order Runge-Kutta with an adaptive step size. In this 
scheme, each integration is performed twice, with two different step 
sizes. For example, if you compute y(x) first using a stepsize of Ax, 
and then using a stepsize of Ax, /2, the difference in the results, €1, 
should be a measure of the local accuracy. Moreover, since the local 
error term on fourth-order Runge-Kutta integration goes as (Az)°, 
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where €o is some desired accuracy, and Azo is the stepsize that 
will yield that accuracy. Thus, to gain additional accuracy from a 





Runge-Kutta method, you can compare your estimate of the local 
accuracy with some desired accuracy. If it is too large, you can try 
again with a smaller stepsize; if it is too small, you can speed things 
up by increasing Az for the next step of the integration. 





